A Multi- Component Phase Field Crystal Model for Structural Transformations in Metal Alloys 
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We present a new phase field crystal model for structural transformations in multi-component alloys. The 
formalism builds upon the two-point correlation kernel developed in Greenwood et al. for describing struc- 
tural transformations in pure materials [Phys. Rev. Lett. 105, 045702 (2010)]. We introduce an effective two- 
point correlation function for multi-component alloys that uses the local species concentrations to interpolate 
between different crystal structures. A simplified version of the model is derived for the particular case of three- 
component (ternary) alloys, and its equilibrium properties are demonstrated. Dynamical equations of motion for 
the density and multiple species concentration fields are derived, and the robustness of the model is illustrated 
with examples of complex microstructure evolution in dendritic solidification and solid-state precipitation. 
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I. INTRODUCTION 

Engineering alloys require the additions of multiple com- 
ponents to achieve desired properties. This, however, makes 
the investigation of their microstructure evolution and defect 
interactions difficult. The properties, and therefore the re- 
sultant behaviour, of alloys can directly be correlated to the 
chemical make up, microstructure and the phase selection 
processes these alloys undergo upon solidification and sub- 
sequent downstream processing, such as thermal treatments. 
In the case of binary alloys, models of solidification pro- 
cesses such as nucleation, free growth and coarsening kinet- 
ics, segregation and second phase formation have been rel- 
atively well developed. However, for multi-component al- 
loys, the complex interactions involved between the different 
chemical species, dislocations and other defects make such 
phenomena far more difficult to study, even with advances in 
characterization techniques such as conventional and high res- 
olution transmission electron microscopy. 

Advances in modelling have significantly improved our un- 
derstanding of the fundamental nature of microstructure and 
phase selection processes. Notable contributions have been 
made using the phase field methodology (PFM), which has 
been successful at examining mesoscale microstructure evo- 
lution over diffusive time scales. The greatest success of the 
PFM has come in the area of solidification |1-7|. The phase 
field concept has gone far beyond its origins. It is now ca- 
pable of describing, through the introduction of various aux- 
iliary fields, a wealth of phenomena such as multiple crystal 
orientations liS nlOJ . multiple components and phases llTTHTSl . 
defect-solute interactions lfT4l . elasticity |[l5][T6l and plastic- 
ity m- 

There has recently emerged an atomic-scale modelling for- 
malism, called the phase field crystal model (PFC) ifTSl [T9l . 
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This method, operates on atomistic length scales and diffu- 
sive time scales and self-consistently incorporates elasticity, 
multiple crystal orientations, grain boundaries, dislocations, 
and the evolution of microstructure on diffusive time scales. 
For both pure materials and binary alloys. Elder and co- 
workers 1 19 1 and Jin and Khachaturyan |20| have shown that 
PFC models can be formally derived from classical density 
functional theory (CDFT), where the order parameter can be 
related to the atomic probability density II2TI . As such, many 
basic microstructure phenomena can be seen as arising self- 
consistently from a simple fundamental theory described by a 
small set of physically motivated parameters. With the ability 
of the PFC density field to also assume disordered states, it is 
also possible to examine amorphous or glassy states Il22ll23l . 
Phase field crystal models are also exceedingly simple to work 
with numerically. The use of coarse graining approaches has 
further shown that PFC-type models can be used as genera- 
tors of traditional phase field models, as well as so-called am- 
plitude models, essentially phase field models with complex 
order parameters. These make it possible to simulate differ- 
ent crystal orientations and defect structures on mesoscopic 
length and time scales Il24lj27l . and also exploit the scaling 
afforded by adaptive mesh refinement 1 28] . 

A weakness of the early PFC models was their inability to 
systematically describe and control complex crystal structures 
and coexistence between them. Greenwood et al. Il29l l30l 
addressed this shortcoming by introducing a class of multi- 
peaked, two-point direct correlation functions that contained 
some of the salient features of CDFT, but retained the simpli- 
fications that gave the original PFC formalism its numerical 
efficiency. This so-called "XPFC" formalism was later ex- 
tended to binary alloys, and applied to phenomena such as 
eutectic solidification and elastic anisotropy |31 1, solute drag 
ll32l . quasi-crystal formation E3]| . solute clustering and pre- 
cipitation mechanisms in simplified Al-Cu alloys |,34J and 3D 
stacking fault structures in FCC crystals 1351 . 

In this paper, we generalize the XPFC formalism of Green- 
wood et. al to the case of A^-component alloys. The approach 
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begins with the truncated CDFT energy functional of an N- 
component system. At the core of our excess free energy are 
the particle interactions of Ref. ||29l [30l . adapted for differ- 
ent structural phases in alloys by making the interaction ker- 
nel a function of the local species concentrations. We com- 
pute the equilibrium properties of our model for the case of a 
ternary alloy and compare the resulting model phase diagram 
to an experimental ternary system. The dynamics of the model 
are then demonstrated in the context of dendritic solidification 
and solid-state precipitation 

The remainder of this paper is organized as follows. We 
begin with deriving the full A^-component XPFC energy func- 
tional in Sectionjnjfrom a simplified, truncated classical den- 
sity function theory of freezing similar to that of Ramakrishan 
and Yussouff |36|. We then derive a second, simplified ver- 
sion of the model that is the iV-component analogue of the 
previous binary XPFC model in literature. Section III calcu- 



lates the equilibrium properties of the model for a particular 
case of a ternary system via isothermal sections of the phase 
diagram. Section IV presents some numerical examples of 



microstructure evolution by simulating dendritic solidification 
and solid-state precipitation. We end with a summary and con- 
clusions. 



II. XPFC ENERGY FUNCTIONALS FOR A-COMPONENT 
ALLOYS 



A general free energy functional for an A^-component al- 
loy is derived starting from the classical density functional 
theory of freezing energy formalism of Ramakrishan and Yus- 
souff 1 36 1, where each alloy component is written in terms of a 
density field pi. The model is re-written in terms of total den- 
sity and concentration variables to make contact with standard 
models used in the description of alloys. The model is then 
collapsed to a simplified form of the free energy, similar to 
the simplified form for the binary XPFC model of Greenwood 
et al. |31|. Finally, equations of motion for the total density 
and each concentration field are presented for both versions of 
the model free energy. 



Deriving an XPFC Energy Functional for A-component 
Systems 



The free energy functional of an A^-component mixture can 
be described by two contributions; a local free energy for 
each of the A^ density fields and an excess free energy due 
to species interactions. The local free energy is treated as an 
ideal energy which drives the density fields to become uni- 
form. The excess contribution drives the density fields to be- 
come periodic by creating minima in the free energy for these 
states. We can write the free energy functional of the mixture 
as 
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where AFid denotes the ideal energy and AF^x is the ex- 
cess energy which accounts for interactions between atoms 
through correlative interactions. This latter term, gives rise to 
structural symmetry, elasticity and interactions between topo- 
logical defects. The constant fc^ is the Boltzmann constant 
and T the temperature. The differential dr = dxdydz. 

The ideal energy, AFi^, gives the entropic contribution for 
an A^-component system. For small density changes from a 
reference density of each component, it is defined as 
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where A^ denotes the number of components, which are de- 
noted as A, B, C,. . . ,etc., pi is the density of component i, and 
p° is the reference density of component i in the liquid phase 
at co-existence. Following previous PFC models lfT9ll . we de- 
fine a total mass density p = pi and the total reference 

mass density as p° = X^f P°- Following Refs. fl9ll27l[3T1. 
we define concentrations as Ci = pi/p and the corresponding 
reference compositions by c" — pH p° ■ Furthermore, for con- 
venience we define a dimensionless mass density of the form 
n — pI p" — \. With these definitions and the conservation 
condition ^ ■ = 1, Eq. (j2j) simplifies to the dimensionless 
form 
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= (n+l)ln(n + l)-?i + (?i + l)2_]ci In 



(3) 



The excess energy takes into account inter-particle interac- 
tions truncated at two-particle, i.e., A-A,B-B,. . . , N-N, A- 
B,. . . ,A-N, ■ ■ ■ interactions. This can be defined as. 
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where C2 represent all combinations of two-particle corre- 
lations, in this work assumed isotropic (i.e., C2 (r, r') — 
C2 (|r — r'D), between the field describing species i and j, 
respectively, where j = A, B,C, . . . , N. We write Eq. (|4]) 
in terms of the reduced density n and compositions c^. As 
in Refs. |fT9ll3TI . we consider only the lowest order contribu- 
tions of the compositions Ci, which vary on length scales much 
larger than the density n, which are periodic on the scale of 
the lattice constant. This allows us to simplify integrals aris- 
ing from Eq. Q, which couple Ci(r') together with n(r') |47|. 
For example. 



dr'C2^dr-r'|)n(r')c,(r') wc,(r) / dr'C2^'(|r-r'|)n(r'). 



To simplify notation, the notation ?i(r') = n' and Ci(r') = 
is used hereafter. With these simplifications and notations, the 
excess energy of Eq.d?]) can be written in terms of the dimen- 
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sionless variables n and {ci} as 



J dr [n a Cj + a Cj — c° cj] J dr' C^J n' 

i,3 



-]-Y, j d-r[nc^+c,-c1] I dv' c;- (5) 

N 

where C'2"' is the Fourier transform of Cj"* (|r — r'|), and satis- 
fies 



C^^(|k|=0)= /dr'C^^(|r-r'| 



(6) 



and where we have introduced the notation C^'' = 
p°C2 {\r — r'\), which is the direct two-point correlation 
function. 

Collecting terms from Eqs. ^ and (jSj gives the complete 
A^-component free energy functional, written in dimension- 
less form. 



p°. The logarithms in the entropy of mixing (Eq. ([8]l) are left 
unexpanded for convenience. Secondly, the terms with cor- 
relation kernels can be simplified by retaining the long wave- 
length behaviour of all compositions c^, where they vary much 
more slowly than n. Following the procedures outlined in 
Refs. Il26ll27l[37l[38l . it can be shown that upon coarse grain- 
ing, all terms containing linear powers of n or n' in Eq. (j?]) 
vanish. Also, terms containing only concentration fields and a 
correlation function give rise to local products of Ci Cj (which 
arise from the k = part of C2 , and look analogous to the 
last term in Eq. (j7]l) and products between their correspond- 
ing gradients. The reader is referred to Appendix|A]for details 
of the coarse graining procedure applied to terms of Eq. (j7]i. 
After some tedious but straightforward algebra, the above ap- 
proximations lead to the following simplified A^-component 
XPFC free energy functional, 

/I n'^ TL"^ Ti^ 
drl y-r;— +x — +^A^^„ix({cJ)(7i+l) 

If 1 ^ 1 

~-n dr'Ce//(|r-r'|)n' + -^«,,Vc..VcJ, (9) 
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The parameters 77, % and uj are constants, the significance of 
which is discussed further below. The Kij are gradient energy 



1 [ dr [71 Ci + Ci — c°] [ dr'C^^c- (7) coefficients associated with compositional interfaces involv 

2 - J J ins Ci and c . For notational convenience, F is used to denot( 
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where Ai^,„ix({ci}) denotes the ideal entropy of mixing. 



N 



AFmix({ci}) = ^ Cj In ^ 



(8) 



Equation (|7]i is the fuU iV-component PEC model in CDFT 
form. When a form for is specified, it can be used di- 
rectly. However, this form is not convenient to make contact 
with other theories and models in the literature. It will be 
transformed into a simpler form in the next section. 



B. Simplified iV-Component XPFC Free Energy 

It is instructive to reduce the model of Eq. (|7]i to a mini- 
mal form that retains the salient features of the original model 
but can also make contact with previous PEC and phase field 
models. To do so, certain simplifications must be made. 

First, an expansion of the ideal free energy term is taken to 
fourth order in the limit of small n, i.e., around the reference 



ing Ci and Cj. For notational convenience, F is used to denote 
^F/kBTp°. 

The parameters 77 and x corresponding to Eq. (|7]i are for- 
mally equal to one, but hereafter will be treated as free pa- 
rameters that can be used to correct the density dependence of 
the ideal free energy away from the reference density p°, i.e., 
to match the bulk free energy to materials properties. Also, it 
was shown in Ref. 1 38 1 that the k = mode of higher-order 
correlation terms in a CDFT expansion will contribute local 
polynomial terms in Ci and n, analogous to the C'2"'(|k| = 0) 
terms of Eq. (j7]l. These terms can be combined with an expan- 
sion of the AFmix term in Eq. (jTji to produce a messy polyno- 
mial expansion of the local free energy in powers of the el- 
ements of {ci} and n. To keep the form of the free energy 
compact, we have found that it is simpler to introduce a pa- 
rameter, Lo, which modifies the mixing free energy from its 
ideal form, away from the reference compositions c°. Finally, 
in the present work, the gradient energy coefficient tensor will 
be assumed to be diagonal for simplicity, i.e., 
i ^ j and na > for all i. 



= for 



The correlation function in Eq. ( 10 1 is too basic to capture 
the properties of very complex alloys -although it can cap- 
ture some properties of simple alloys. Guided by the form 
of the first term on the second line of Eq. (jTjl, it can be seen 
that higher-order correlation functions will contribute terms of 



the form Ci Cj Ck C^J^ , Ci Cj Ck ciC/''"'', etc. To emulate such 
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higher-order non-local contributions effectively, we introduce 
an effective correlation function of the form 



where the following shorthand notations have been made. 
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The Xi are as yet undetermined polynomial functions of the 
elements of {cj}. The role of the Xi is to determine the re- 
sultant local crystalline structure by interpolating between the 
kernels C2 (defined below), which define the base equilib- 
rium crystal structures of each pure component i. The interpo- 
lation is done through appropriately constructed polynomial 
expansions of the elements of {cj}. The order of Xi depends 
on the number of components in the system and can be made 
as high as required to smoothly interpolate from one correla- 
tion kernel to another. We have found that Eq. (Ill, through 
appropriate choices of Xi, combined with other model param- 
eters, is robust enough to model a wide variety of alloy sys- 
tems. 

The model in Eq. (|9]l captures the usual features of other 
PFC models, while allowing for a very easy control of a wide 
range of crystal structures in different phases. It is motivated 
from considerations of classical density functional theory but 
simplified enough to make numerically tractable simulations 
possible, as will be shown below. Finally, we note that the 



form of the expansion in Eq. ( 11 1 is dimensionally motivated 



from higher-order terms in CDFT but is flexible enough to 
model experimentally relevant multi-component alloys quan- 
titatively using, for example, thermodynamic databases. 



C. Dynamics 



Equations of motion for the density n and each of the con- 
centration fields Cj follow conserved dissipative dynamics. 
Namely the dimensionless density n obeys 



(12) 



+ c«, 

while the dynamics of each composition field, Ci, evolve ac- 
cording to 
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The coefficients Af„ and Mc- denote the mobility of the den- 
sity and each concentration, respectively, and strictly speak- 
ing can be functions of the fields. The noise terms C,n 
and model coarse grained thermal fluctuations on den- 
sity and concentrations Ci, respectively. They formally sat- 
isfy (C,(r,i)Cg(r',<')> = -AV'xaiv - v')6{t - t'), where 
q denotes the density or one of the concentration fields, A (x 
MqksT and Xa(r — r') is the inverse Fourier transform of 
a Gaussian function, which, following Tegze and co-workers 
Ref. If39l . can be generalized to have a high frequency cut off 
for frequencies above 27r / a, where a is the lattice constant. 
The precise form of A, which sets the scale of the thermal 
fluctuations is not properly understood in the context of PFC 
modelling but is the object of several investigations. In the 
applications illustrated in this paper, the noise is left out of 
simulations. 



III. TERNARY SYSTEMS 

In this section, we reduce the simplified free energy func- 
tional of section (IIBi to the case of three-components, or 
ternary alloys. We first describe the ternary free energy func- 
tional, followed by a discussion of the effective correlation 



function chosen for ternary systems. With the free energy and 
effective correlation in hand, we demonstrate the equilibrium 
properties of our model by calculating the ternary phase dia- 
grams for a generic A-B-C system and a simplified Al-Cu- 
Mg system. 



A. Free Energy Functional 

Speciahzing Eq. (|9]l for 3-components, denoted here as A, 
B and C, reduces it to 

= j dvl--7^-+x^+ u^^F::1 (n + 1) (15) 
-^nj dr'C:}j{\r-r'\)n' + f\VcA 



where 
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and 
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+ Xb{ca, cb) ^(|r - r'l) + Xc{ca. cb) C^"" {\v 



(17) 
r'l). 



In arriving at Eq. ( 15 i, the conditions cc = ^ — ca ^ cb and 



c°Q = 1 — — Cg have been used, and the cross gradient 
concentration terms in A and B have been neglected. 

The effective ternary correlation kernel, Cfjj, is defined by 
Xi such that Xa + Xb + Xc = 1 at all compositions, anal- 
ogous to the case for the XPFC binary model 1,31 1 . Their par- 
ticular form, is chosen here to model the generic properties 
of eutectic systems. However, by careful alteration of other 
parameters, other alloy systems can be modelled, e.g. isomor- 
phous and peritectic systems ifSTl . Here the Xi used are, 

Xa{ca,cb) = 'Sc\ + "^cacb - 2c^ - 2c\cb - '2cac% 
Xb{ca,cb) = 2cacb + 3c^B - '^c\cb - 2cac% ~ 2c| 



XcicA,CB) = 1 - 3c 



24 
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2c% - ACACB 
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bility of equilibrium structures is typically unchanged f30l. 
The situation is similar for alloys, where the k = mode will 
have a negligible contribution for phases that remain relatively 
close to the reference density. Deviations of phases away from 
the reference density will be manifested in the average density 
dimension of the phase diagram. Here, it is assumed that the 
average density rio = to simplify the demonstration of the 
model. Of course, the more complex situations where both the 
concentration and average density need to be modelled can be 
treated by adding suitable k = contributions, or by choos- 
ing the appropriate coefficients in the bulk free energy. Thus, 
without loss of generality, we will assume no additional con- 
stant to the correlation function here. 



C. Ternary Dynamics 



For the case of 3-component alloys, the dynamical equa- 



tions of motions in Eqs.( 12 1-( 13 1 reduce to 
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mix 
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B. Correlation Functions C2' 

The XPFC formalism is best suited for numerical simula- 
tion in Fourier space. The pure component correlation func- 
tions C2'(|r — r'l) are thus constructed directly in Fourier 
space, where they are denoted 6*2* (k). Each component, i, 
contributes a correlation function that supports the desired 
equilibrium crystal structure for a pure component. A Fourier 
space peak of C'2'(k) |30|, for a given mode, j, is denoted by 



c: 
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The total correlation function for component i, CI 



by the envelope of all peaks C 



2r 



is defined 
The first exponential in 



Eq. ( 19l sets the temperature scale via a Debye-Waller pref- 
actor that employs an effective temperature parameter, a. We 
also define an effective transition temperature, (TMj, which 
subsumes the effect of planar and atomic densities associated 
with the family of planes corresponding to mode j |31 1. The 
second exponential sets the position of the reciprocal space 
peak at kj, which defines the inverse of the interplanar spacing 
for the j'^ family of planes in the equilibrium unit cell struc- 
ture of component i. Each peak is represented by a Gaussian 
function, with aj being the width of the peak, j. The {aj} 
have been shown in Ref. Il30l to set the elastic and surface 
energies, as well as their anisotropic properties. 

It is noted that the k = mode of all correlation functions 
is essentially zero. In principle, as discussed above, the k = 
mode of these correlation functions can have their effects im- 
plicitly reflected through local coefficients in the free energy. 
In the case of a pure material, a nonzero peak height at k = 
in the correlation function merely shifts the local free energy 
at densities away from the reference density, however the sta- 



dcA 
dt 



1)^^ - T^n—^ n-AC^V^CAk 

OCA 2 dCA I 



where Af„, Mc^ and Mc^ are dimensionless mobility coeffi- 
cients for density and compositions fields. They are set to 1 
here, since it is the intent of this paper to introduce the model 
and its physical features. 



D. Equilibrium Properties 

Ternary equilibrium is defined by co-existence of bulk 
phases, e.g. solids-solid^, liquid-solid^-solid^, etc. The gov- 
erning properties, e.g. partitioning, of such an equilibrium 
state can be determined from standard thermodynamic min- 
imization methods. In general, for 3-component alloys, free 
energy minimization is defined by a common plane tangent to 
the free energy wells of any two or three coexisting phases. 
This construction is a geometrical representation of the state- 
ment that the chemical potentials and grand potentials of any 
two phases are equal with respect to each component. Here, 
we construct isothermal ternary phase diagrams by examin- 
ing all combinations of phase coexistence (e.g., solido, -liquid, 
solids-solid^, etc.). Procedures for calculating phase dia- 
grams for PFC models are well-documented ifTSl [T9l |29l |3TI 
and the approach used here will only be summarized. 

For solid phases, the density field, which varies on atomic 
length scales, is approximated using a multi-mode approxi- 
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mation given by 



1. Generic Ternary Eutectic Alloy 



1=1 



27r 



(21) 



where is the lattice spacing of the solid phase i and Ni de- 
notes the number of mode famihes (families of planes) in the 
unit cell of phase z, Aj is the amplitude associated with the 
j"^ family of planes. Each mode contains Nj reciprocal lat- 
tice peaks, enumerated by the index I. Strictly speaking, there 
is a distinct amplitude, Ai^j, for each reciprocal lattice peak. 
However, for the purposes of simplifying the construction of 
phase diagrams (i.e., working with the fewest number of vari- 
ables to minimize), they are assumed constant leading to Aj. 
Each index I in the family j has a corresponding reciprocal 
lattice vector k/ j, normalized by the lattice spacing. 



Substituting Eq. (21 1 into Eq. (15 1, and integrating over one 
unit cell, the free energy can be calculated for each phase as a 
function of ca, cb and the amplitudes Aj. Since amplitudes 
are non-conserved fields, the resulting free energy is then min- 
imized with respect to each Aj |31 1. The result is substituted 
back into the free energy. After this procedure, we are left 
with a free energy landscape J^^'^^;(cyi, c_b), where J^^^i rep- 
resents an amplitude-minimized solid free energy. In keep- 
ing with the discussion of the previous sections, we assume 
that the average density of all phases is close to the reference 
density, i.e., Uo = 0. For the liquid phase, the free energy 
J^^^(c^, cs) is trivially computed by setting all Aj — 0\ 



With the free energy landscapes of liquid and solids, the 
phase boundary lines between a combination of phases at a 
given temperature parameter, cr, are computed by solving the 
following set of equations simultaneously. 



^^'cA 



A first example of the equilibrium properties of the ternary 
XPFC model are demonstrated for a system where all three 
components {A, B, and C) are structurally similar, differ- 
ing only in their equilibrium lattice spacings. Here two- 
dimensional (2D) square symmetry is assumed as the equi- 
librium structure for each pure component, which in this con- 
text implies that all C2 have the same number of peaks, with 
the corresponding ratios of their positions in reciprocal space 
being the same. However, each structure is differentiated by 
the absolute positions (kj) of each peak. Though it has not 
been done in this initial work, by adjusting the widths (aj) of 
each peak, each element can also be differentiated by different 
elastic and surface energies. The full list of parameters used 
to construct the phase diagrams in this subsection are listed in 
the caption of Fig. [T] 

Allowing all three components to have square structural 
symmetry, at sufficiently low temperature we can construct 
a bulk solid free energy landscape describing multiple solid 
phases, described by an effective lattice parameter (a^") that 
is a weighted average of the individual lattice parameters 
of all three components, using the interpolation functions of 
Eq.( 18 I, namely, a'*"" — qaXa + 03X3 + ac'Xc. This 



(22) 



leads to the solid-liquid free energy landscape in Fig.[TJa) for 
a = 0.17, where the values of all other parameters are speci- 
fied in the figure caption. The corresponding isothermal phase 
diagram is illustrated in Fig. [ijc), which is constructed form 
the coexistence lines calculated between the liquid phase and 
the different solid-solution phases, using the set of conditions 
in Eq. (22 1. Figure [TJb) shows an isothermal cut at a higher 
temperature, i.e., a — 0.182, depicting an increased region 
where the bulk liquid is stable compared to the solid phases. 
At sufficiently low temperature, the free energy admits eutec- 
tic coexistence of three phases. We construct an isothermal cut 
right above the eutectic temperature, i.e., — 0.164, shown 
in Fig.[T|d). The corresponding concentrations ca and cb in 
Fig. [T] are given as fractions, where unity represents pure A 
or B, respectively, along each axis of the phase diagram. 



2. Simplified Al-Cu-Mg Type Alloy 



where the last of these implies. 



i^i^A 



(23) 



The superscripts / and J denote any two phases in equi- 
librium (e.g. liquid-solida), respectively. The expressions 
/i^^ ~ dj^ /dcA and ^l.^ — df^ /Ocb are the chemical po- 
tentials of phase / with respect to the concentrations ca and 
Cb, respectively, with analogous definitions for ji'l^ and fi"^^. 
The expressions ft^ and il'^ are the grand potentials of phases 
/ and J, respectively. See Appendix [B] for further details on 



calculating phase diagrams. The set of conditions in Eq. (22 1 



along with Eq. (B3i defining the average concentration, can 
be solved to find the four equilibrium concentrations (two per 
phase) defining coexistence on a given tie line. 



The parameters of the ternary XPFC model can be chosen 
to produce sections of experimental phase diagrams qualita- 
tively, as in the work of Fallah et al. |40|, where the present 
ternary model is used to model precipitation in a 2D represen- 
tation of the Al-Cu-Mg system. Here, we demonstrate how the 
equilibrium properties of a portion of the Al-rich (simplified) 
part of the Al-Cu-Mg phase diagram can be described quanti- 
tatively by the ternary XPFC model. An experimental phase 
diagram at 400°C is shown in Fig.|2|b), taken from Ref. |41 1. 

Consider the part of the phase diagram for (Al)— /S-^ out- 
lined by the red dashed line and circled solid phases in the 
experimental phase diagram shown in Fig. |2|b), and ignor- 
ing the (AI)h-S and (A1)h-/3h-T phase regions. In the dilute-Mg 
region, a eutectic transition occurs between the Al-rich, (Al)- 
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FIG. 1: (Colour online) Ternary eutectic system: (a) Solid and liquid energy landscapes of a square-square-square (A-B-C) system at temper- 
ature parameter a = 0.17. Corresponding phase diagrams at temperatures (b) a = 0.182, (c) cr = 0.17 and (d) a — 0.164. The parameters for 
ideal free energy and entropy of mixing were rj = 1.4, x ~ ^ ~ 0.005, while reference concentrations were c'X ~ 0.333 and c% = 0.333. 
Widths of the correlations peaks are taken an = 0.8 and aio — V^an for all phases (required for isotropic elastic constants in a solid phase 
with square symmetry Bil l. The peak positions for the given structures are knA ~ (81/38)7r and kioA ~ V^kuA for a, kns = (54/29)7r 
and fcios = v^^^iis for P and fcnc = Stt and fcioc = V^^iic for 7. The effective transition temperatures are set to a^jj — 0.55 for all 
family of planes in all phases. The concentrations on the isothermal phase diagrams are read in a Cartesian coordinate system. 



fee phase, and an intermediate phase which has a tetrag- 
onal crystal structure. The eutectic system of (AVj-d has a 
small solubility for Mg, however past the maximum solubil- 
ity limit, there exists other intermediate phases terminating 
at the cubic /3-phase. The equilibrium lattice constants (and 
thus the positions of the reciprocal space peaks) of 9 and (3 
phases are determined by interpolating between those of Al 
with 32.5 at.%Cu, and Al with 38.5 at.%Mg, respectively. For 
simplicity, we assume a square structural symmetry for all 
three equilibrium phases, and like the preceding section, the 
effective lattice constant is interpolated by weighting by local 
solute compositions, ccu and cmq- The parameters (77, x, u), 
along with the peak widths aj are chosen to give a satisfactory 
mapping of the solubility limits of the (Al)-phase for Cu and 
Mg to those in the experimental phase diagram of Fig. |2|b), 
for a range of temperature parameters (cr). The full list of 
parameters used to construct the phase diagrams in this sub- 
section are listed in the caption of Fig. [2] 

Figure |2|a), shows the free energy landscape for the solid 
at cr = 0.04. Figure [2|c) shows the coiTesponding isother- 
mal phase diagram at cr = 0.04, where the inset shows a 
zoomed in image of the Al-rich corner Comparing the inset 



with the experimental phase diagram, reasonable agreement 
is evident between the calculated and the experimental phase 
diagram sections. Figure [2];d) shows the isothermal phase di- 
agram for a = 0.155. At this higher temperature (still below 
the eutectic), there is an increase in the solubility limits of the 
phase boundaries. Section [TV B| will use this phase diagram to 
demonstrate solid-state precipitation. 



IV. APPLICATIONS 

The binary XPFC approach was previously demonstrated 
as a tool with which to model the role of defects and elasticity 
in structural phase transformations that operate over diffusive 
time scales. Further to these capabilities, the ability to have 
multi-component interactions between solute atoms and de- 
fects now makes it possible to examine much more complex 
interactions of the above atomic-scale effects with different 
solutes, and their diffusion. This capability opens a myriad 
of possibilities for applications for microstructure engineer- 
ing in materials. This section showcases some applications 
of the XPFC multi-component model presented in this work. 
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FIG. 2: (Colour online) Al-Cu-Mg Phase diagram: (a) Solid and liquid energy landscapes of a square-square-square ((Al)— /3 — 6) system at 
temperature a = 0.04, (b) The Al-rich side of an isothermal cut (at 400°C) from the experimental phase diagram of the Al-Cu-Mg system 
taken from Ref. |41|. Dashed circles mark the regions of the Al-rich (Al), Cu-rich {6} and Mg-rich (/3) regions considered for reconstruction 
by the model phase diagram. Reconstructed phase diagrams at temperatures (c) a = 0.04 and (d) a — 0.155. The parameters for ideal free 
energy and entropy of mixing were ry = 1.4, x = 1^ = 0.005, Cq^ — 0.333 and c'^fg = 0.333. Widths of the correlations peaks are 
ail ~ O.Sandaio — \/2qii for all phases. The peak positions are fcii(Ai) = 2tv, fcio(A!) = V^kn^^Ai), kng = (2.0822)7r, fcioe — \/2kiig, 
fcii/9 ~ (1.8765)7r and kiop = \/2A;ii/3. For all family of planes, aMj ~ 0.55, in all phases. The maxima in concentrations ccu and CMg 
are rescaled from unity, 1, to correspond to the Cu and Mg-content in the ^-phase and /3-phase given by the experimental phase diagram, i.e., 
~ 32.5 and ~ 38.5 at.%, respectively. The concentrations on the isothermal phase diagrams are read in a Cartesian coordinate system. 



In particular, using the phase diagrams from the previous sec- 
tion, we demonstrate dendritic solidification and precipitation 
in the presence of ternary components. These phenomena are 
paradigms of microstructure evolution of relevance to mate- 
rials engineering applications and are strongly influenced by 
diffusion of impurities, elastic strain, crystal anisotropy and 
defect structures. 



A. Dendritic Solidiflcation 

Dendritic solidification arises when a supercooled liquid is 
quenched into the solid-liquid coexistence part of the phase 



diagram. Figure ([3]) shows snapshots in time of a dendritic 
crystal in a ternary alloy. The simulation was done using the 
phase diagram in Fig. [TJb). Simulations were conducted in a 
2D domain of size 768a x 768a, where a is the lattice spac- 
ing. A uniform grid spacing and discrete time of Aa; — 0.125 
(which makes for a domain of size 6144 x 6144 grid points) 
and At = 3 were used and equations of motion, Eq. (20i, 



were solved semi-implicitly in Fourier space. The initial con- 
ditions consisted of a small circular seed of diameter, d — 8a 
of 7-phase, seeded in liquid at a temperature of ct = 0.182. 
The initial concentration of solute components A and B was 
uniform in both phases and set to the values ca — 0.1 and 
cb = 0.1. Several time slices of the simulation domain, show- 
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FIG. 3; (Colour online) Early-time dendritic solidification in a 
ternary alloy, simulated using the phase diagram of Fig. [T|b). The 
quench temperature is cr = 0.182 and the initial solute compositions 
are uniform and set to the alloy averages, ca = 0.1 and cb = 0.1. 
Each Column of images represents a different time during the sim- 
ulation. The times shown are: (a) 1000 (b) 3000 and (c)7000 itera- 
tions. From bottom to top, each row displays the progression of n, 
Cb and ca, respectively, with ca plotted in the colour range from 
white (lowest concentration) to red (highest concentration) and cb is 
plotted in the colour range from white (lowest concentration) to blue 
(highest concentration). 

ing the fields (n, ca and cb) at early tines, are shown in Fig. [3] 
As time progresses during the simulation. Fig. [3|a)-(c), 
dendritic growth is evident. The crystal develops a character- 
istic 4-fold symmetry of the underlying square crystal struc- 
ture, produced with the coiTelation function for the given pure 
component of the 7-phase. The top two rows show the time 
evolution of the concentration fields (from left to right), ca 
and Cb, respectively, indicating the interface boundary layer 
for each component. Both solutes, A and B, reach their max- 
imum solute content at the interface of the growing dendrite, 
in agreement with the solute rejection mechanism of crystal 
growth. The bottom row shows the evolution of the density. 
There is also evidence of the associated density jump at the 
interface between solid and liquid phases as depicted by the 
light halo like region around the interface. Figure |4] shows 
a composite view of the dendrite at later time, highlighting 
in each quadrant one of the three fields. This simulation de- 
picts multiple diffusing species, density changes and surface 
tension anisotropy. In a larger numerical domain (where mul- 
tiple dendrites can be grown), grain boundaries would also 
naturally emerge. It is noteworthy that these physical ingre- 
dients arise self-consistently and are very straightforward to 
simulate numerically. We also note that side-branching of the 
growing dendrite is not observed in Fig. |3]due to the size of 
the simulation domain and the exclusion of thermal noise in 




FIG. 4: (Colour online) Dendritic solidification at time, t = 10000, 
displaying the values of all three fields. The top left and bottom right 
quadrants show the density field, n. The inset on the bottom right, is 
a zoomed-in image of the rectangular area marked in black, revealing 
the structure of the underlying square lattice and density wave struc- 
ture through the interface. Top right quadrant shows ca, with colour 
range white (low concentration) to red (high concentration), while 
bottom left shows cb, with colour range white (low concentration) 
to red (high concentration). Both show the high solute content at the 
interface. 

the dynamical equations. 



B. Solute Clustering and Precipitation 

Many properties of engineering alloys are typically at- 
tained through downstream processing following solidifica- 
tion. These downstream processes typically involve either 
thermo-mechanical manipulations or heat treatment of the as- 
cast microstructure. One of the most important aims is to in- 
duce certain phase transformations in the as-cast primary solid 
matrix to help strengthen alloys, a process known as precipi- 
tation hardening. In this subsection we demonstrate this pro- 
cess using the ternary XPFC model developed in this work. 
In particular, we illustrate the initial stages of a heat treatment 
process leading to solute clustering, the precursor stage of pre- 
cipitation in Al-Cu-Mg alloys. The details of this process have 
been reported elsewhere 140J. 

Solute clustering/early-stage precipitation simulations were 
preformed using the equilibrium properties calculated for the 
(Al)— /3-6' system in Figs. |2|c) and |2|d). Simulations were 
peiformed on a 2D rectangular mesh with giid spacing Aa; = 
0.125 and time step At = 10. Dynamical equations were 
solved semi-implicitly in Fourier space. Initial conditions 
consisted of distorted single-phase structures, through the in- 
troduction of a uniform distribution of dislocations, and a uni- 
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form composition everywhere of ccu = 11 and CMg = 0.2 
at.%. All simulations were initially solutionized for some time 
at (T = 0.155, following which they were quenched/aged at a 
temperature a — 0.04. During ageing, small clusters initially 
appear with higher Mg and/or Cu-content than that of the ma- 
trix. As time progresses, some of these clusters decrease in 
size and Cu-Mg-content, or vanish entirely. A few, however, 
stabilize, as shown by the typical stabilized clusters "a" and 
"6" in Fig. |5]:a)-(c) and (d)-(f) for Al-l.lCu and Al-l.lCu- 
0.2Mg alloys, respectively. In contrast, for either alloy, when 
we increase the ageing temperature within the single-phase 
(Al) region , e.g., a = 0.145, no clustering is observed and 
the initial distortions are removed from the matrix. 

Experiments in quenched/aged Al-Cu and Al-Cu-Mg al- 
loys l,42i - i44J have found increasing evidence that the inter- 
action of ternary impurities and quenched-in defects such as 
dislocations [45] dynamically reduce the local nucleation bar- 
rier for precipitation at locations in the matrix. We have also 
found that the addition of Mg into an Al- 1 . 1 Cu alloy promotes 
clustering and refinement of the final microstructure, as seen 
in the simulation data of Fig. (j5]l. The clustering phenomenon 
observed in these simulations can be attributed to the propen- 
sity for solute segregation to defects and surrounding areas to 
relieve stresses induced by the presence of said defects, in this 
case dislocations. As more solute aggregates to dislocations, 
the size of the cluster increases but the structural nature of the 



cluster also begins to approach that of the next nearest stable 
solid phase. As this process continues and the ever growing 
cluster attracts more solute, it creates additional stresses in the 
surrounding matrix. This in turn draws nearby dislocations 
to the cluster in attempts to relieve these additional stresses 
caused by solute accumulation. An extensive investigation 
of solute clustering mechanisms, in presence of quenched- 
in bulk crystal defects, has been done through a quantitative 
analysis of the system energetics in binary alloys in ll34l and 
recently in ternary alloys, using the present model POl . 



V. SUMMARY 

This paper reported a new phase field crystal model for 
structural phase transformations (XPFC) in multi-component 
alloys. The details of the model derivation were discussed. 
A simplified version of the model was specialized for ternary 
alloys and its equilibrium properties were shown. The dy- 
namics of the model were demonstrated on two phenomena 
of relevance to microstructure evolution in materials science. 

This is the first multi-component PFC model, and as such is 
able to capture the complex kinetics of solidification and elas- 
tic and plastic effects on solid state processes, such as cluster- 
ing and precipitate growth. This model has been used in a sep- 
arate work 1 40 1 to support recent experiments on the elusive 
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mechanisms of the early stages of clustering and precipitation. 

The phase field crystal methodology was introduced to 
create a bridge between the atomic and traditional phase 
field regimes. As a relatively novel method, many works 
in this area of materials science are working to validate the 
physics of PFC models. As the first phase field crystal 
model for A^-component alloys, this work has demonstrated 
some important thermodynamic and kinetic properties of the 
model. Moreover, aside from the model's quantitative and 
self-consistent nature, it is particularly simple to operate nu- 
merically. It is expected that this model can thus be used to 
elucidate the role of multiple solutes in phenomena governed 
by atomic-scale elasticity and defects operating on diffusional 
time scales. 
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Appendix A: Long Wavelength Limit 



where we define 



c,(k)^ / dr'c,(r') 



(A4) 



Considering the long wavelength limit, we take a Taylor series 
expansion of the correlation function in powers of around 
k = 0. This results in. 



Q = 



(A5) 



i,j 1=0 



5(k2)' 



k r 



Cj(k)e"" 



J 

k=0 



We note that to invoke the long wavelength limit, we could 
have also Taylor expanded the concentration, Cj (r'), at r' = r 
as is done in Refs. Il25l - l27l or employed the multi-scale expan- 
sion used in Refs. 1371 |38l . All these methods, though differ- 
ent and require different mathematical treatments, are found 
to be equivalent. Retaining, to lowest order, terms up to order 
Z = 1, we have 



1 ^ 



kr 



k=0 



9(k2 



cj(k)e'''-'". (A6) 



k=0 



In section II B we reduced the free energy in Eq. (j7]i into 
a simplified multi-component XPFC energy functional. In the 
process of doing this, we simplified terms by considering the 
long wavelength limit where the concentration varies much 
more slowly than the density field. This Appendix details the 
steps of how some terms of Eq. (j7]i can be simplified to derive 
the simplified free energy functional in Eq. (|9|. 



1. Terms Coupling Product of ct and Cj with C2 

We begin first with terms involving a coupling of two con- 
centration fields with a correlation function. As a concrete 
example, consider the term 

TV 

Q = -\T.j drc,{r) I dr'C^^"(|r-r'|)c,(r'), (Al) 



in Eq. (|7]l, where we have used the more explicit notation for 
clarity. (The other terms follow analogously.) To proceed, we 
rewrite the correlation function in a Fourier series of the form. 



(A2) 



Substituting Eq. ( |A2| l into Eq. ( |A1[ ) yields. 



N 

Q = -\Y.j drc^r) I rfkC'^^"(|k|)c,(k)e-''-, (A3) 



Using the definition of the inverse Fourier transform, we 



recast Eq. (A6 1 as 



G ^-^Yllij I drc,{r) Cj{r) 

1 ^ r 



where we have used the following definitions. 



and 



k=0 



a(k2) 



(A8) 



(A9) 



k=0 



It is thus clear that the first term in Eq. ( A7 1 will contribute 
terms that renormalize the coefficient of the terms in the 
entropy of mixing, if Eq. (|9jl were expanded about Cj = c°. 
In this work, the 7^^ terms are neglected, and their role is sub- 
sumed in an effective manner, for convenience, into the pref- 
actor uj in Eq. (|9|. The second term in Eq. ( A7 1 can be recast 



into gradient energy terms analogous to those used in Cahn- 
Hilliard or Ginzburg-Landau theories. To do so, we perform 
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integration by parts, yielding, 



where are the reciprocal lattice vectors and Am(r) are 
slowly varying amplitudes corresponding to each reciprocal 



Q = -\Y1^^o j d'^c^{r) Cj(r) 

+ 2 ^ «»j / dr Vc,(r) • Vcj(r). (AlO) 
Lastly, for clarity of exposition, we separate the gradient terms 



n 



in Eq. ( AlOi, which yields 



lattice vector, m. Substituting Eq. ( A14i into Eq. (A13 1 gives 

N 

i,j m 

N 

2 E^».^ E / *<c,(r)V2 (A„(r)e''i'" '-) . 

(A15) 



1 ^ /■ 

+ 2E'*" y |Vc,(r)|' 

+ 2 E / Vc,(r) • Vcj(r). 



Expanding the Laplacian in Eq. ( A15 1 gives. 



(All) 



N 

^ - -2E^'^:'E / drc°c,{r)A„,{r)e'^-- 

i,j m 

1 ^ f 

-oE^^jE / rfrc°c,(r)e'i'"'-/:™^„(r), (A16) 



where 



2iqm • V — is a covariant opera- 



The second term of Eq. ( |A11| | gives gradients terms of the tor that assures rotational invariance of the free energy in the 



Cahn-Hilliard form, while the third line yields cross terms. In 
this work, for simplicity, we set kij = for i 7^ j. We note 
that such cross terms can become important when studying 
certain phenomena and/or when higher-order alloying inter- 
actions are considered. 



2. Correlation Kernels Containing Linear Terms in n 

To demonstrate the long wavelength limit of terms linear in 
density in Eq. (|7]i, we consider, as an example, the term 

N 

'H = ^\Y. Idrc'ic,{v) I dv'Cl^{\v-v'\)n{v'). (A12) 
Substituting the Fourier series expansion of the correlation 



long wavelength limit. It is noted that each term in Eq. ( A16 1 
only contains one rapidly oscillating variable, i.e., e"i'" ''. If 
we apply the so-called "quick and dirty" |,46J analogue of the 
volume averaging method employed in Refs. Il26l IZTl (which 
amounts to decoupling slowly varying fields inside integrals 
from rapidly varying phase factors, thus making the integrals 
effectively vanish when integrated over one unit), we obtain 

It is straightforward to show that all other terms in Eq. (|7]) 
that are linear in n, such as, 

TV 

n^-lY^j dvn{v)c,{v) I dr'C^^(|r-r'|)c,(r'), 

id 

(A17) 

similarly vanish upon coarse graining. It should also be evi- 



function, Taylor expanding the correlation as in Eq. ( A5 1 (re- 
taining the lowest order term), and taking the inverse Fourier 
transform yields. 



1 ^ 

?i = --^7ij / drc°Cj{r)n{r) 



1 ^ r 

+ 2 E '^'J J dvc1c,{v) (-V2) n(r), (A13) 

i,3 



dent from the above considerations that if Eq. ( A12 1 contained 
an ri(r) • • • n(r') combination, then Eq. (AI61 would contain 



terms with phase factors of different combinations of sums of 
two reciprocal lattice vectors. Some of these two-vector com- 
binations would add up to zero causing their corresponding 
terms to survive upon integration. 



3. Volume Averaging 



where 7^^ and are defined by Eqs. ( A8 1 and ( A9 1, respec- 
tively. 



Equation ( A16 1 can more formally be analyzed using a vol- 



The density, n(r), in Eq. ( A13 1 is rapidly varying. Its lead- 
ing order representation is defined by a single-mode approxi- 
mation of the form 



(A14) 



ume averaging convolution operator |26|, defined by 

(/(r))y = ^-/ dr7(r')xy(r-r'), (A18) 

where /(r') is the function being course grained and V is the 
coarse graining volume. The function xv in the integrand of 



Eq. ( A18 1 is a smoothing function that is normalized to unity. 
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I.e., 



tween two phases. However for multi-component systems, 
for the present case of a ternary alloy (represented by solute 



/OO ^ ^ ^ i w 

drxvi^ — r') — 1 (A19) compositions A and B), the set of conditions in Eq. (22i are 

-OO under-determined and cannot uniquely define all phase con- 



A convenient form of xv is given by 

1 (r-r'):^ 

Xv{r-r') = ^-e . 



(A20) 



under-determined and cannot uniquely define all phase con- 
centrations. This is because, for a ternary system at a speci- 
fied temperature and pressure, there is not generally a single 
tie line which specifies phase boundaries between coexisting 
phases but, rather multiple tie lines defining the boundary be- 
tween any two phases. 



In the long wavelength limit, L ^ a where L ^ V^^'^, 

in d-dimensions, while Lc is the length scale of variation of 
the concentration field. This condition implies that the func- 
tion Xv (r) varies on dimensions much larger than the lattice 
constant a — 27r/|qm| but much less then the length scale of 
variation of the concentration, Ci{r). Equation ( A18 i defines 
a noninvertible limiting procedure that can be used to average 
a function over some volume. 

It is instructive to apply the volume averaging procedure 
to the first term in Eq. ( |A16 



The under-determined set of conditions in Eq. ( 22 1 contain 



variables c^, 



-A ' '-B 



and Cg in phases / and J respectively. 



To close this system, an additional condition is necessary to 
provide a fourth equation relating the concentrations. A con- 
venient fourth condition is the lever rule, which relates weight 
fractions of phases to the average concentration. For clarity, 
we specify it here for ternary solid (a) and liquid (L) phases. 



For convenience we define 
(p[Y) = c°Cj(r)A,„(r). It is noted that (/)(r) varies on scales 
much larger than the lattice constant since it is comprised of 
slowly varying functions. Using the definition of 0(r), the 



and 



cb 



(Bl) 



(B2) 



first integral of Eq. ( A16 1 can be written as 



where ca and cb are the average alloy compositions for com- 
ponents A and B respectively and a;^ and Xa represent the 
equilibrium volume fractions of liquid and a respectively, and 
satisfy x^ + Xa = 1. Combining this last relation between 
the volume fractions and Eqs. ( |B1| | and ( |B2[ ) gives, the last 
equihbrium condition. 



t.j m 



CA 



CB 



(A21) 



(B3) 



Since (/)(r') varies more slowly than the scale of variation of 
Xv, it is reasonable to expand it in a Taylor series about r' = 
r. Substituting (f>{r') ~ 0(r) — V(/>(r) • (r — r') into the above 
expression leads to 



Equation ( 22 1 together with Eq. ( B3 1 comprises a complete set 
of equations which can admit unique tie line solutions, i.e., 
solutions for c^,c'^,Cg and eg in the solid-liquid example 
just considered. 

With the free energy functions generally being highly non- 
linear, it is not possible to find analytical solutions to Eqs. (|22l) 



N f / f Imear, It IS not possible to find analytical solutions to Eqs. (|22|) 

?iv=--E7>jE / drf /dr'xylr-rOe''!" ""' and (|B3]l, and they must be solved numerically. One approach 

V \*i to <5nppifv thp tpmnpratiirp and thpn ra*;tpr thrniiah thp nha-sp 



W(j){r)- / dr'xy(r-r') (r - r') e' 



(A22) 



The noninvertible procedure was introduced in the second line 
of Eq. ( A21 1. In the lon g wav elength limit, when |qm|i — > 
OO, both integrals in Eq. ( A22 1 vanish as ^ (|qmL|) ^, mak- 



is to specify the temperature and then raster through the phase 
space of average concentrations ca and cb, where the ras- 
tering is done by taking discrete steps in steps of A ca and 
Acs, respectively (for practical purposes its convenient to 
set AcA — Acb = Ac). For each pair of ca and cb. 



ing Hv similarly vanish. 



Appendix B: Phase Diagram Calculation 



Equation ( 22 1 provides a system of equations that are ex- 



act when one needs to determine the equilibrium properties 
of a given system. For a binary system, where the number of 



equations in Eq. ( 22 1 is reduced by one, at a specified tem- 
perature and pressure they are sufficient to specify exactly the 
unique phase concentrations corresponding to the tie line be- 



Eqs. (22 1 and (B3 1 can be solved numerically. The solutions 
yield c'^, Cg and c%. A unique solution for each pair of ca 
and Cb defines one tie line. The collection of all such tie lines 
maps out the coexistence phase boundaries between any two 
phases, in the case considered here, L and a. Where no so- 
lutions are admitted correspond to single phase regions where 
no tie lines exist. It is expected that the smoothness of the 
phase boundaries, when plotted for graphical purposes, will 
depend on the step size, A c, chosen to discretize the average 
concentration values. 

The above mentioned recipe can still require intensive com- 
putation, requiring a solution of four equations in four un- 
knowns for combinations of average concentration pairs 
(where M is the discretized number of average concentration 
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values for a given component). Since this paper is intended 
to demonstrate the main features of our new multi-component 
(demonstrated for a ternary) PFC model, we adopted a sim- 



pler approach to compute the phase diagrams in section 111 D 



In particular, we fixed one of the equilibrium concentrations 
in Eq. (|22]), assuming it is a valid solution at that temperature. 
We then solved for the remaining three unknown concentra- 



tions using Eq. ( [22) l, repeating this M times, once for each 
discrete value of the selected equilibrium concentration. Fixed 
concentrations were rastered in steps of A c. Once again, a 
unique solution defines a tie line between coexisting phases, 
say L and a. If no solutions exist, we are in single phase 
regions where no tie lines exist. 
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